home *** CD-ROM | disk | FTP | other *** search
/ PC Electronics Plus 3 / PC Electronics Plus 3.iso / coil01 / mathl.c < prev    next >
C/C++ Source or Header  |  1992-12-20  |  7KB  |  366 lines
  1. extern double MACHEP, MAXNUM;
  2.  
  3. /* Factorials of integers from 0 through 33 */
  4. double fac[] = {
  5.   1.0,
  6.   1.0,
  7.   2.0,
  8.   6.0,
  9.   2.4E1,
  10.   1.2E2,
  11.   7.2E2,
  12.   5.04E3,
  13.   4.032E4,
  14.   3.6288E5,
  15.   3.6288E6,
  16.   3.99168E7,
  17.   4.790016E8,
  18.   6.2270208E9,
  19.   8.71782912E10,
  20.   1.307674368E12,
  21.   2.0922789888E13,
  22.   3.55687428096E14,
  23.   6.402373705728E15,
  24.   1.21645100408832E17,
  25.   2.43290200817664E18,
  26.   5.109094217170944E19,
  27.   1.12400072777760768E21,
  28.   2.585201673888497664E22,
  29.   6.2044840173323943936E23,
  30.   1.5511210043330985984E25,
  31.   4.03291461126605635584E26,
  32.   1.0888869450418352160768E28,
  33.   3.04888344611713860501504E29,
  34.   8.841761993739701954543616E30,
  35.   2.6525285981219105863630848E32,
  36.   8.22283865417792281772556288E33,
  37.   2.6313083693369353016721801216E35,
  38.   8.68331761881188649551819440128E36
  39. };
  40.  
  41.  
  42. /*                            ellpe.c
  43.  *
  44.  *    Complete elliptic integral of the second kind
  45.  *
  46.  *
  47.  *
  48.  * SYNOPSIS:
  49.  *
  50.  * double m1, y, ellpe();
  51.  *
  52.  * y = ellpe( m1 );
  53.  *
  54.  *
  55.  *
  56.  * DESCRIPTION:
  57.  *
  58.  * Approximates the integral
  59.  *
  60.  *
  61.  *            pi/2
  62.  *             -
  63.  *            | |                 2
  64.  * E(m)  =    |    sqrt( 1 - m sin t ) dt
  65.  *          | |    
  66.  *           -
  67.  *            0
  68.  *
  69.  * Where m = 1 - m1, using the approximation
  70.  *
  71.  *      P(x)  -  x log x Q(x).
  72.  *
  73.  * Though there are no singularities, the argument m1 is used
  74.  * rather than m for compatibility with ellpk().
  75.  *
  76.  * E(1) = 1; E(0) = pi/2.
  77.  *
  78.  *
  79.  * ACCURACY:
  80.  *
  81.  *                      Relative error:
  82.  * arithmetic   domain     # trials      peak         rms
  83.  *    DEC        0, 1       13000       3.1e-17     9.4e-18
  84.  *    IEEE       0, 1       10000       2.1e-16     7.3e-17
  85.  *
  86.  *
  87.  * ERROR MESSAGES:
  88.  *
  89.  *   message         condition      value returned
  90.  * ellpe domain      x<0, x>1            0.0
  91.  *
  92.  */
  93.  
  94. /*
  95. Cephes Math Library, Release 2.1:  February, 1989
  96. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  97. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  98. */
  99.  
  100. static double PE[] = {
  101.   1.53552577301013293365E-4,
  102.   2.50888492163602060990E-3,
  103.   8.68786816565889628429E-3,
  104.   1.07350949056076193403E-2,
  105.   7.77395492516787092951E-3,
  106.   7.58395289413514708519E-3,
  107.   1.15688436810574127319E-2,
  108.   2.18317996015557253103E-2,
  109.   5.68051945617860553470E-2,
  110.   4.43147180560990850618E-1,
  111.   1.00000000000000000299E0
  112. };
  113. static double QE[] = {
  114.   3.27954898576485872656E-5,
  115.   1.00962792679356715133E-3,
  116.   6.50609489976927491433E-3,
  117.   1.68862163993311317300E-2,
  118.   2.61769742454493659583E-2,
  119.   3.34833904888224918614E-2,
  120.   4.27180926518931511717E-2,
  121.   5.85936634471101055642E-2,
  122.   9.37499997197644278445E-2,
  123.   2.49999999999888314361E-1
  124. };
  125.  
  126. double ellpe(x)
  127. double x;
  128. {
  129. double polevl(), log();
  130.  
  131.  
  132. if( (x <= 0.0) || (x > 1.0) )
  133.     {
  134.     if( x == 0.0 )
  135.         return( 1.0 );
  136.     printf( "ellpe domain error\n" );
  137.     return( 0.0 );
  138.     }
  139.  
  140. return( polevl(x,PE,10) - log(x) * (x * polevl(x,QE,9)) );
  141. }
  142.  
  143. /*                            ellpk.c
  144.  *
  145.  *    Complete elliptic integral of the first kind
  146.  *
  147.  *
  148.  *
  149.  * SYNOPSIS:
  150.  *
  151.  * double m1, y, ellpk();
  152.  *
  153.  * y = ellpk( m1 );
  154.  *
  155.  *
  156.  *
  157.  * DESCRIPTION:
  158.  *
  159.  * Approximates the integral
  160.  *
  161.  *
  162.  *
  163.  *            pi/2
  164.  *             -
  165.  *            | |
  166.  *            |           dt
  167.  * K(m)  =    |    ------------------
  168.  *            |                   2
  169.  *          | |    sqrt( 1 - m sin t )
  170.  *           -
  171.  *            0
  172.  *
  173.  * where m = 1 - m1, using the approximation
  174.  *
  175.  *     P(x)  -  log x Q(x).
  176.  *
  177.  * The argument m1 is used rather than m so that the logarithmic
  178.  * singularity at m = 1 will be shifted to the origin; this
  179.  * preserves maximum accuracy.
  180.  *
  181.  * K(0) = pi/2.
  182.  *
  183.  * ACCURACY:
  184.  *
  185.  *                      Relative error:
  186.  * arithmetic   domain     # trials      peak         rms
  187.  *    DEC        0,1        16000       3.5e-17     1.1e-17
  188.  *    IEEE       0,1        30000       2.5e-16     6.8e-17
  189.  *
  190.  * ERROR MESSAGES:
  191.  *
  192.  *   message         condition      value returned
  193.  * ellpk domain       x<0, x>1           0.0
  194.  *
  195.  */
  196.  
  197.  
  198. /*
  199. Cephes Math Library, Release 2.0:  April, 1987
  200. Copyright 1984, 1987 by Stephen L. Moshier
  201. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  202. */
  203.  
  204. static double PK[] =
  205. {
  206.  1.37982864606273237150E-4,
  207.  2.28025724005875567385E-3,
  208.  7.97404013220415179367E-3,
  209.  9.85821379021226008714E-3,
  210.  6.87489687449949877925E-3,
  211.  6.18901033637687613229E-3,
  212.  8.79078273952743772254E-3,
  213.  1.49380448916805252718E-2,
  214.  3.08851465246711995998E-2,
  215.  9.65735902811690126535E-2,
  216.  1.38629436111989062502E0
  217. };
  218.  
  219. static double QK[] =
  220. {
  221.  2.94078955048598507511E-5,
  222.  9.14184723865917226571E-4,
  223.  5.94058303753167793257E-3,
  224.  1.54850516649762399335E-2,
  225.  2.39089602715924892727E-2,
  226.  3.01204715227604046988E-2,
  227.  3.73774314173823228969E-2,
  228.  4.88280347570998239232E-2,
  229.  7.03124996963957469739E-2,
  230.  1.24999999999870820058E-1,
  231.  4.99999999999999999821E-1
  232. };
  233. static double C1 = 1.3862943611198906188E0; /* log(4) */
  234.  
  235.  
  236.  
  237. double ellpk(x)
  238. double x;
  239. {
  240. double polevl(), log();
  241.  
  242.  
  243.  
  244. if( (x < 0.0) || (x > 1.0) )
  245.     {
  246.     printf( "ellpk domain error\n" );
  247.     return( 0.0 );
  248.     }
  249.  
  250. if( x > MACHEP )
  251.     {
  252.     return( polevl(x,PK,10) - log(x) * polevl(x,QK,10) );
  253.     }
  254. else
  255.     {
  256.     if( x == 0.0 )
  257.         {
  258.         printf( "ellpk singularity\n" );
  259.         return( MAXNUM );
  260.         }
  261.     else
  262.         {
  263.         return( C1 - 0.5 * log(x) );
  264.         }
  265.     }
  266. }
  267.  
  268.  
  269. /*                            polevl.c
  270.  *                            p1evl.c
  271.  *
  272.  *    Evaluate polynomial
  273.  *
  274.  *
  275.  *
  276.  * SYNOPSIS:
  277.  *
  278.  * int N;
  279.  * double x, y, coef[N+1], polevl[];
  280.  *
  281.  * y = polevl( x, coef, N );
  282.  *
  283.  *
  284.  *
  285.  * DESCRIPTION:
  286.  *
  287.  * Evaluates polynomial of degree N:
  288.  *
  289.  *                     2          N
  290.  * y  =  C  + C x + C x  +...+ C x
  291.  *        0    1     2          N
  292.  *
  293.  * Coefficients are stored in reverse order:
  294.  *
  295.  * coef[0] = C  , ..., coef[N] = C  .
  296.  *            N                   0
  297.  *
  298.  *  The function p1evl() assumes that coef[N] = 1.0 and is
  299.  * omitted from the array.  Its calling arguments are
  300.  * otherwise the same as polevl().
  301.  *
  302.  *
  303.  * SPEED:
  304.  *
  305.  * In the interest of speed, there are no checks for out
  306.  * of bounds arithmetic.  This routine is used by most of
  307.  * the functions in the library.  Depending on available
  308.  * equipment features, the user may wish to rewrite the
  309.  * program in microcode or assembly language.
  310.  *
  311.  */
  312.  
  313.  
  314. /*
  315. Cephes Math Library Release 2.1:  December, 1988
  316. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  317. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  318. */
  319.  
  320.  
  321. double polevl( x, coef, N )
  322. double x;
  323. double coef[];
  324. int N;
  325. {
  326. double ans;
  327. int i;
  328. double *p;
  329.  
  330. p = coef;
  331. ans = *p++;
  332. i = N;
  333.  
  334. do
  335.     ans = ans * x  +  *p++;
  336. while( --i );
  337.  
  338. return( ans );
  339. }
  340.  
  341. /*                            p1evl()    */
  342. /*                                          N
  343.  * Evaluate polynomial when coefficient of x  is 1.0.
  344.  * Otherwise same as polevl.
  345.  */
  346.  
  347. double p1evl( x, coef, N )
  348. double x;
  349. double coef[];
  350. int N;
  351. {
  352. double ans;
  353. double *p;
  354. int i;
  355.  
  356. p = coef;
  357. ans = x + *p++;
  358. i = N-1;
  359.  
  360. do
  361.     ans = ans * x  + *p++;
  362. while( --i );
  363.  
  364. return( ans );
  365. }
  366.